Aimant permanent

2.a. Charges magnétiques

Dans cette partie, on considère le problème général du champ magnétique créé par un aimant permanent, en l'absence de tout courant électrique.

Les matériaux ferromagnétiques peuvent posséder une aimantation en l'absence de source de champ extérieur. Un aimant permanent est défini comme un volume (V) en tout point duquel la densité volumique de moment magnétique $\overrightarrow{M}$ (l'aimantation), est supposée connue et constante.

Le champ d'aimantation est équivalent à un vecteur densité de courant électrique égal au rotationnel de l'aimantation ([1]). En régime quasi stationnaire, le champ magnétique créé par l'aimant obéit donc à l'équation de Maxwell-Ampère à laquelle ce courant équivalent est ajouté :

$$\overrightarrow{rot}\overrightarrow{B}={\mu }_0(\overrightarrow{J}+\overrightarrow{rot}\overrightarrow{M})\text{(1)}$$

où $\overrightarrow{J}$ est la densité de courant électrique, qui est nulle dans le problème considéré ici.

L'équation $\mathrm{(1)}$ est une version de l'équation de Maxwell-Ampère (en régime quasi stationnaire) valable dans la matière et faisant intervenir des champs moyens, c'est-à-dire des champ moyennés dans l'espace sur une échelle beaucoup plus grande que l'échelle atomique

La propriété de conservation du flux magnétique reste valable pour le champ moyen :

$$div\overrightarrow{B}=0\text{(2)}$$

On définit le vecteur excitation magnétique par :

$$\overrightarrow{H}=\frac{\overrightarrow{B}}{{\mu }_0}-\overrightarrow{M}\text{(3)}$$

L'équation $\mathrm{(1)}$ s'écrit alors :

$$\overrightarrow{rot}\overrightarrow{H}=\overrightarrow{J}\text{(4)}$$

et la conservation du flux magnétique :

$$div\overrightarrow{H}=-div\overrightarrow{M}\text{(5)}$$

Le vecteur excitation magnétique $\overrightarrow{H}$ s'exprime en A/m. Il est parfois nommé champ magnétique et dans ce cas le vecteur $\overrightarrow{B}$ est désigné comme le vecteur densité de flux magnétique, qui est une dénomination tout à fait correcte, ou bien induction magnétique, qui est plus discutable. Dans les problèmes de magnétostatique dans le vide, le vecteur $\overrightarrow{H}$ n'a pas d'intérêt et il est d'usage de désigner le vecteur $\overrightarrow{B}$ comme étant le champ magnétique. Nous préférons garder cette dénomination et désigner $\overrightarrow{H}$ par excitation magnétique.

En l'absence de courant électrique, les deux équations vérifiées par l'excitation magnétique sont :

$$\begin{array}{cc}\hfill & \overrightarrow{rot}\overrightarrow{H}=\overrightarrow{0}\hfill \\ \hfill & div\overrightarrow{H}=-div\overrightarrow{M}\hfill \end{array}$$

Ces deux équations sont similaires à celles vérifiées par le champ électrostatique. Cela nous amène à définir une densité de charge magnétique (par analogie avec la densité volumique de charge électrique) par :

$${\rho }_m=-div\overrightarrow{M}\text{(6)}$$

Plus précisément, si le champ électrostatique est analogue à l'excitation magnétique, ρ_m et analogue à ρ/ε₀.

Soit S la surface fermée qui délimite le volume V de l'aimant. En principe, l'aimantation, qui est un champ moyen défini à l'échelle mésoscopique, devrait varier continûment lorsqu'on passe de l'intérieure de l'aimant à l'extérieur. Cependant, la distance sur laquelle se fait la variation de l'aimantation au voisinage de la surface est négligeable à l'échelle macroscopique. Il peut donc être intéressant, pour des raisons de simplicité du modèle, de considérer que l'aimantation est non nulle dans le volume V et nulle en dehors de ce volume (dans l'air qui entoure l'aimant). En conséquence, il subit une discontinuité sur la surface S et sa divergence n'est donc pas définie sur cette surface. Ce problème est résolu en considérant l'intégrale de la densité de charge magnétique sur le volume :

$${\iiint }_V{\rho }_{m}dv=-{\iiint }_{V}div\overrightarrow{M}dv=-{\iint }_S\overrightarrow{M}\cdot \overrightarrow{n}dS\text{(7)}$$

où $\overrightarrow{n}$ est le vecteur normal sortant du volume V. Une densité surfacique de charge magnétique doit donc être introduite sur la surface de l'aimant, définie par :

$${\sigma }_m=\overrightarrow{M}\cdot \overrightarrow{n}\text{(8)}$$

où $\overrightarrow{n}$ désigne la normale de la surface de l'aimant dans le sens sortant.

Par analogie avec le champ électrostatique, nous pouvons écrire l'excitation magnétique en un point P de l'espace sous la forme d'une intégrale sur le volume de l'aimant plus une intégrale sur sa surface :

$$\overrightarrow{H}\left(P\right)=\frac{1}{4\pi }{\iiint }_V\frac{\overrightarrow{QP}}{\left|\right|\overrightarrow{QP}|{|}^3}{\rho }_{m}dv+\frac{1}{4\pi }{\iint }_S\frac{\overrightarrow{QP}}{\left|\right|\overrightarrow{QP}|{|}^3}{\sigma }_{m}dS\text{(9)}$$

où Q désigne un point quelconque du volume ou de la surface. Cette relation permet de calculer l'excitation magnétique générée par un ou plusieurs aimants permanents. Cependant, elle ne permet pas de la calculer en présence de parties ferromagnétiques qui ne sont pas aimantées de manière permanente (par exemple les noyaux de bobines en fer doux) car l'aimantation de ces parties n'est pas connue a priori.

Pour un barreau aimanté, on considère généralement que l'aimantation est uniforme dans le barreau. On supposera donc que l'aimantation est uniforme dans l'aimant. Dans ce cas, il reste seulement l'intégrale sur la surface :

$$\overrightarrow{H}\left(P\right)=\frac{1}{4\pi }{\iint }_S\frac{\overrightarrow{QP}}{\left|\right|\overrightarrow{QP}|{|}^3}\overrightarrow{M}\cdot \overrightarrow{n}dS\text{(10)}$$

L'excitation en un point P étant connue, le champ magnétique s'en déduit par :

$$\begin{array}{ccc}\hfill & \overrightarrow{B}={\mu }_0(\overrightarrow{H}+\overrightarrow{M})dansl\text{'}aimant\hfill & \hfill \text{(11)}\\ \hfill & \overrightarrow{B}={\mu }_0\overrightarrow{H}dansl\text{'}air\hfill & \hfill \text{(12)}\end{array}$$

L'équation $\mathrm{(10)}$ permet de faire un calcul direct de l'excitation magnétique par un calcul d'intégrale de surface.

Une autre manière d'aborder le problème consiste à introduire un potentiel magnétique. En effet, le rotationnel de l'excitation magnétique étant nul en l'absence de courants électriques, il existe un potentiel ${\Phi }_m$ tel que :

$$\overrightarrow{H}=-\overrightarrow{grad}{\Phi }_m\text{(13)}$$

L'opérateur divergence appliquée à cette équation conduit à l'équation :

$$\Delta {\Phi }_m=-{\rho }_{m}dansl\text{'}aimant,\Delta {\Phi }_m=0dansl\text{'}air\text{(14)}$$

où Δ est l'opérateur laplacien. Il s'agit de l'équation de Poisson, similaire à l'équation de Poisson pour le potentiel électrostatique. Cette équation n'a de sens que si la charge magnétique est définie, autrement dit si la divergence de l'aimantation est définie en tout point. Il faudra donc adopter un modèle qui respecte cette condition, ce qui exclue l'introduction d'une charge magnétique surfacique.

Il faut aussi écrire les conditions limites à la frontière entre l'air et l'aimant. L'équation $div\overrightarrow{B}=0$ implique la continuité de la composante normale du champ magnétique. Dans le cas où l'aimantation est discontinue à l'interface aimant-air, si $\overrightarrow{n}$ désigne le vecteur normal dirigé de l'aimant vers l'air, cette condition s'écrit :

$$({\overrightarrow{H}}^{aimant}+\overrightarrow{M})\cdot \overrightarrow{n}={\overrightarrow{H}}^{air}\cdot \overrightarrow{n}\text{(15)}$$

La composante normale de l'excitation magnétique a donc une discontinuité égale à l'aimantation sur la surface.

Si l'aimantation varie continûment entre l'aimant et l'air, la continuité de la composante normale du champ magnétique et de l'aimantation implique la continuité de la composante normale de l'excitation magnétique.

L'équation $\overrightarrow{rot}\overrightarrow{H}=\overrightarrow{0}$ implique la continuité de la composante tangentielle :

$$\overrightarrow{n}\wedge ({\overrightarrow{H}}^{aimant}-{\overrightarrow{H}}^{air})=\overrightarrow{0}\text{(16)}$$

L'équation de Poisson peut être résolue numériquement avec une condition limite de potentiel nul sur les bords du domaine, qui doit être de très grande taille par rapport à l'aimant. Une méthode de discrétisation par différence finie associée à une méthode de résolution multigrille devrait constituer un moyen très efficace de résoudre numériquement cette équation.

2.b. Barreau aimanté cylindrique

Le barreau aimanté est cylindrique, d'axe (Oz). Sa section droite peut être circulaire ou rectangulaire. On suppose que l'aimantation dans le volume du barreau est uniforme et dirigée selon l'axe (Oz) :

$$\overrightarrow{M}=M{\overrightarrow{u}}_z\text{(17)}$$

On s'intéresse aux faces sud et nord du barreau, respectivement situées en z=z_s et z=z_n car c'est dans leur voisinage qu'une charge magnétique est présente.

3.a. Hypothèses

On s'intéresse au cas particulier d'un aimant ou de deux aimants cylindriques de section circulaire, alignés sur le même axe.

Figure pleine page

Un point P est repéré par ses coordonnées cylindriques (r,θ,z). Tout plan contenant l'axe (Oz) est un plan de symétrie de l'aimantation (vecteur axial). En conséquence, c'est aussi un plan de symétrie de l'excitation magnétique (vecteur axial). Il y a par ailleurs invariance par rotation autour de l'axe. Il s'en suit que :

$$\begin{array}{ccc}\hfill & \overrightarrow{H}=H_r(r,z){\overrightarrow{u}}_r+H_z(r,z){\overrightarrow{u}}_z\hfill & \hfill \text{(24)}\\ \hfill & {\Phi }_m(r,z)\hfill & \hfill \text{(25)}\end{array}$$

3.b. Calcul du champ magnétique

Le calcul direct de l'excitation magnétique consiste, pour chaque face des aimants, à calculer l'intégrale $\mathrm{(18)}$. Il faut calculer l'intégrale pour chaque face (2 faces pour un aimant, 4 faces pour deux aimants) et sommer les champs obtenus pour obtenir le champ complet. Voyons comme se fait le calcul de l'intégrale pour une face d'un aimant de rayon a, c'est-à-dire un disque de rayon a. La densité de charge magnétique sur cette face est σ_m=M ou σ_m=-M selon qu'il s'agit de la face nord ou de la face sud de l'aimant d'aimantation M.

L'intégrale s'écrit :

$$\overrightarrow{H}\left(P\right)=\frac{{\sigma }_m}{4\pi }{\iint }_{disque}\frac{\overrightarrow{QP}}{\left|\right|\overrightarrow{QP}|{|}^3}dS\text{(26)}$$

Le rapport $\overrightarrow{H}/M$ est sans dimensions. En conséquence, le rayon et la longueur de l'aimant peuvent être exprimés dans l'unité que l'on souhaite, par exemple le millimètre.

Pour expliciter cette intégrale, on peut considérer que θ=0. Notons (r₁,θ₁,z₁) les coordonnées cylindriques du point Q situé sur une face.

On a :

$$\overrightarrow{QP}=(r-r_{1}cos\left({\theta }_1\right)){\overrightarrow{u}}_x-r_{1}sin\left({\theta }_1\right){\overrightarrow{u}}_y+(z-z_1){\overrightarrow{u}}_z\text{(27)}$$

L'élément de surface s'écrit dS=r₁dθ₁dr₁. Les deux composantes de l'excitation magnétique s'écrivent donc :

$$\begin{array}{ccc}\hfill & H_r(r,z)=\frac{{\sigma }_m}{2\pi }{\int }_0^{a}d{r}_1{\int }_0^{\pi }d{\theta }_1\frac{r_1(r-r_{1}cos\left({\theta }_1\right))}{(r^2+r_1^2-2r{r}_{1}cos\left({\theta }_1\right)+{(z-z_1)}^2{)}^{\frac{3}{2}}}\hfill & \hfill \text{(28)}\\ \hfill & H_z(r,z)=\frac{{\sigma }_m}{2\pi }{\int }_0^{a}d{r}_1{\int }_0^{\pi }d{\theta }_1\frac{r_1(z-z_1)}{(r^2+r_1^2-2r{r}_{1}cos\left({\theta }_1\right)+{(z-z_1)}^2{)}^{\frac{3}{2}}}\hfill & \hfill \text{(29)}\end{array}$$

Pour calculer numériquement ces intégrales, on doit effectuer un maillage du demi-disque de rayon a. L'angle θ₁ est échantillonné de la manière uniforme :

$$\begin{array}{ccc}\hfill & \Delta \theta =\frac{\pi }{N_{\theta }}\hfill & \hfill \text{(30)}\\ \hfill & {\theta }_{1,i}=i\Delta \theta i=0,1,\dots N_{\theta }-1\hfill & \hfill \text{(31)}\end{array}$$

Pour le rayon r₁, il est préférable de choisir un échantillonnage qui permette d'avoir des mailles d'aire constante entre le centre et les bords du disque. Celui-ci devrait convenir :

$$\begin{array}{ccc}\hfill & \Delta r=\frac{a}{\sqrt{N_r}}\hfill & \hfill \text{(32)}\\ \hfill & r_{1,j}=\sqrt{j}\Delta rj=0,1,\dots N_r-1\hfill & \hfill \text{(33)}\end{array}$$

La figure suivante représente la maille (i,i+1,j,j+1) :

Figure pleine page

Son aire est :

$$\Delta S=\frac{\Delta \theta }{2}(r_{1,j+1}^2-r_{1,j}^2)=\frac{\Delta \theta {\left(\Delta r\right)}^2}{2}(j+1-j)=\frac{\Delta \theta {\left(\Delta r\right)}^2}{2}\text{(34)}$$

Toutes les mailles ont bien la même aire.

Les approximations des intégrales sont obtenues en considérant que l'intégrale sur une maille est égale à l'aire de la maille multipliée par la valeur de la fonction intégrée au centre de la maille.

$$H_r(r,z)\approx \frac{{\sigma }_m}{2\pi }\frac{\Delta \theta {\left(\Delta r\right)}^2}{2}\sum _{i=0}^{N_{\theta }-1}\sum _{j=0}^{N_r}\frac{(r-r_{1,j+\frac{1}{2}}cos\left({\theta }_{1,i+\frac{1}{2}}\right))}{(r^2+r_{1,j+\frac{1}{2}}^2-2r{r}_{1,j+\frac{1}{2}}cos\left({\theta }_{1,j+\frac{1}{2}}\right)+{(z-z_1)}^2{)}^{\frac{3}{2}}}\text{(35)}$$ $$H_z(0,z)\approx \frac{{\sigma }_m}{2\pi }\frac{\Delta \theta {\left(\Delta r\right)}^2}{2}\sum _{i=0}^{N_{\theta }-1}\sum _{j=0}^{N_r}\frac{(z-z_1)}{(r^2+r_{1,j+\frac{1}{2}}^2-2r{r}_{1,j+\frac{1}{2}}cos\left({\theta }_{1,j+\frac{1}{2}}\right)+{(z-z_1)}^2{)}^{\frac{3}{2}}}\text{(36)}$$

Il peut être intéressant de valider ce calcul en comparant H_z à son expression sur l'axe (Oz), qui admet une forme analytique :

$$H_z(0,z)=\frac{2\pi {\sigma }_m}{4\pi }(z-z_1){\int }_0^a\frac{r_1}{(r_1^2+{(z-z_1)}^2{)}^{3/2}}d{r}_1=\frac{{\sigma }_m}{2}[\frac{z-z_1}{|z-z_1|}-\frac{z-z_1}{(a^2+{(z-z_1)}^2{)}^{1/2}}]\text{(37)}$$

De part et d'autre de la surface, on a :

$$\begin{array}{ccc}\hfill & H_z(0,z_1^{+})=\frac{\sigma }{2}=\frac{M}{2}\hfill & \hfill \text{(38)}\\ \hfill & H_z(0,z_1^{-})=-\frac{\sigma }{2}=-\frac{M}{2}\hfill & \hfill \text{(39)}\end{array}$$

ce qui confirme la discontinuité de la composante normale de l'excitation magnétique, égale à l'aimantation M.

Numpy permet de faire des calculs de somme double très rapidement. Pour cela, on utilise tout d'abord la fonction numpy.meshgrid, qui permet de crééer une grille à partir de deux tableaux. Par exemple, pour créer une grille définie par x_i=i et y_j=j avec 3 points sur chaque axe, on écrit :

import numpy as np
x = [0,1,2]
y = [0,1,2]
xx,yy = np.meshgrid(x,y)
                   

print(xx)
--> array([[0, 1, 2],
       [0, 1, 2],
       [0, 1, 2]])

print(yy)
--> array([[0, 0, 0],
       [1, 1, 1],
       [2, 2, 2]])

Le tableau xx contient les valeurs de x pour les points de la grille, le tableau yy contient les valeurs de y. L'évaluation d'une fonction a deux variables sur cette grille se fait alors de la manière suivante (pour une fonction universelle) :

def f(x,y):
    return x+y
z = f(xx,yy)
                   

print(z)
--> array([[0, 1, 2],
       [1, 2, 3],
       [2, 3, 4]])

La somme des valeurs sur les nœuds de la grille se calcule aisément :

print(z.sum())
--> 18

Les fonctions suivantes effectuent le calcul des sommes $\mathrm{(35)}$ et $\mathrm{(36)}$ :

def HrHz(r,z,z1,rr,tt,deltaS,M):
    b2 = (z-z1)**2
    cos = np.cos(tt);
    D = (r**2+rr**2-2*r*rr*cos+b2)**1.5
    Hr = (r-rr*cos)/D
    Hz = (z-z1)/D
    Hr = Hr.sum()*deltaS/(2*np.pi)*M
    Hz = Hz.sum()*deltaS/(2*np.pi)*M
    return (Hr,Hz)
    
def grille(Nr,Ntheta,a):
    delta_theta = np.pi/Ntheta
    delta_r = a/np.sqrt(Nr)
    theta1 = (0.5+np.arange(Ntheta))*delta_theta
    r1 = np.sqrt(0.5+np.arange(Nr))*delta_r
    deltaS = delta_theta*delta_r**2*0.5
    rr,tt = np.meshgrid(r1,theta1)
    return (rr,tt,deltaS)
               
def champH(Nr,Ntheta,r,z,z1,a,M):
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    return HrHz(r,z,z1,rr,tt,deltaS,M)
                 

Voici un exemple :

Nr = 100
Ntheta = 50
a = 1
z = 1
z1 = 0
r = 0
M = 1
(Hr,Hz) = champH(Nr,Ntheta,r,z,z1,a,M)

                 

print((Hr,Hz))
--> (0.0, 0.14644532315842412)

Pour comparaison, la fonction suivante calcule le champ sur l'axe, défini par $\mathrm{(37)}$ :

def Hz_axe(z,z1,a,M):
    return M/2*((z-z1)/abs(z-z1)-(z-z1)/(a**2+(z-z1)**2)**0.5)

Hz = Hz_axe(z,z1,a,M)
                 

print(Hz)
--> 0.14644660940672627

Pour faire une comparaison plus complète, traçons H_z en fonction de z pour r=0 :

from matplotlib.pyplot import *
z = np.linspace(0.01,5,100)
Nr = 100
Ntheta = 50
a = 1
z1 = 0
r = 0
M = 1
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
Hz1 = np.zeros(len(z))
for k in range(len(z)):
    Hr,Hz1[k] = HrHz(r,z[k],z1,rr,tt,deltaS,M)
Hz2 = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz1,"r-",label='approché')
plot(z,Hz2,"k",label='exact')
grid()
xlabel('z')
ylabel('Hz')
legend(loc='upper right')
ylim(0,1)
                 

Hz-axe-1face.pdf

Pour cette taille de grille, le calcul approché de l'intégrale donne un résultat convenable sauf à petite distance du disque. À très petite distance du disque, la fonction à intégrer prend des valeurs très grandes lorsque le point Q se trouve proche de P et elle varie beaucoup lorsqu'il s'éloigne de cette position. Pour améliorer la précision du calcul de l'intégrale, il faut augmenter la finesse de la subdivision de r :

z = np.logspace(-4,-2,10)
Nr = 1000
Ntheta = 50
a = 1
z1 = 0
r = 0
M = 1
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
Hz1 = np.zeros(len(z))
for k in range(len(z)):
    Hr,Hz1[k] = HrHz(r,z[k],z1,rr,tt,deltaS,M)
Hz2 = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz1,"r-",label='approché')
plot(z,Hz2,"k",label='exact')
grid()
xlabel('z')
ylabel('Hz')
legend(loc='upper right')
xscale('log')
                 

Hz-axe-1face-2.pdf

Une augmentation d'un facteur 100 du nombre de points selon r permet bien d'augmenter la précision, mais à très courte distance, il y a toujours une divergence de l'intégrale.

Il faut donc adapter le nombre N_r à la distance z et ne pas faire le calcul à une distance trop petite. Il faut trouver un critère pour choisir ce nombre. On peut par exemple considérer que le rayon de la maille centrale doit être de l'ordre de z/p, où p est un paramètre à ajuster. Cette condition conduit à :

$$N_r=E\left[\right(\frac{pa}{\left|z\right|}{)}^2]\text{(40)}$$

mais N_r ne doit pas descendre en dessous d'une certaine valeur.

z = np.linspace(0.01,5,100)
Ntheta = 50
a = 1
z1 = 0
r = 0
M = 1
Hz1 = np.zeros(len(z))
p = 3
Nr_min = 10
for k in range(len(z)):
    Nr = max(int((p*a/abs(z[k]))**2),Nr_min)
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    Hr,Hz1[k] = HrHz(r,z[k],z1,rr,tt,deltaS,M)
Hz2 = Hz_axe(z,z1,a,M)
figure()
plot(z,Hz1,"r-",label='approché')
plot(z,Hz2,"k",label='exact')
grid()
xlabel('z')
ylabel('Hz')
legend(loc='upper right')
ylim(0,1)
                 

Hz-axe-1face-3.pdf

Il reste à savoir si la méthode fonctionne aussi lorsque le point P est proche du disque mais hors de l'axe. L'aire de la maille qui se trouve au plus proche du point P a la même aire quelle que soit la position de P. On peut donc raisonnablement s'attendre à un calcul approché correct également dans ce cas.

Traçons les composantes H_r et H_z ainsi que la norme de H, en fonction de r, pour une distance z petite :

r = np.linspace(0,2,50)
z = 1e-1
a = 1
z1 = 0
M = 1
Hz1 = np.zeros(len(r))
Hr1 = np.zeros(len(r))
Nr_min = 10
p = 2
Nr = max(int((p*a/abs(z))**2),Nr_min)
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
for k in range(len(r)):    
    Hr1[k],Hz1[k] = HrHz(r[k],z,z1,rr,tt,deltaS,M)
H1 = np.sqrt(Hr1**2+Hz1**2)
Hz2 = np.zeros(len(r))
Hr2 = np.zeros(len(r))
p = 3
Nr = max(int((p*a/abs(z))**2),Nr_min)
(rr,tt,deltaS) = grille(Nr,Ntheta,a)
for k in range(len(r)):    
    Hr2[k],Hz2[k] = HrHz(r[k],z,z1,rr,tt,deltaS,M)
H2 = np.sqrt(Hr2**2+Hz2**2)
figure()
subplot(311)
plot(r,Hr1,'r-')
plot(r,Hr2,'b-')
ylabel('Hr')
grid()
subplot(312)
plot(r,Hz1,'r-')
plot(r,Hz2,'b-')
ylabel('Hz')
grid()
subplot(313)
plot(r,H1,'r')
plot(r,H2,'b-')
ylabel('H')
xlabel('r')
grid()
                  

HrHz-prox-1face.pdf

Comme souvent en calcul numérique, le moyen le plus simple de vérifier que le calcul est assez précis est de faire varier le paramètre qui contrôle la précision, ici p, et de comparer les résultats. Les deux courbes précédentes sont calculées pour p=2 et p=3. La valeur p=2 semble suffisante.

Dans l'objectif de tracer des lignes de champ, la méthode de calcul directe semble intéressante, à condition de démarrer les lignes de champ à une distance pas trop petite des faces des aimants. À grande distance, le calcul précis se fait avec une faible complexité.

Le calcul du champ créé par un aimant est fait par la fonction suivante, qui renvoie les composantes de $\overrightarrow{H}$ et de $\overrightarrow{B}$. L'aimant est caractérisé par son rayon a, les côtes de ses faces z_s et z_n et son aimantation M.

def champAimant(r,z,a,zs,zn,M,mu0=1):
    p = 3
    Nr_min = 10
    Ntheta = 50
    Nr = Nr_min
    if abs(r)<=a:
        Nr = max(int((p*a/abs(z-zs))**2),Nr_min)
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    (Hr1,Hz1) = HrHz(r,z,zs,rr,tt,deltaS,-M)
    Nr = Nr_min
    if abs(r)<=a:
        Nr = max(int((p*a/abs(z-zn))**2),Nr_min)
    (rr,tt,deltaS) = grille(Nr,Ntheta,a)
    (Hr2,Hz2) = HrHz(r,z,zn,rr,tt,deltaS,M)
    Br1 = mu0*Hr1
    Br2 = mu0*Hr2
    if abs(r)<a and z>min(zs,zn) and z<max(zs,zn):
        Bz1 = mu0*(Hz1+M)
        Bz2 = mu0*(Hz2+M)
    else:
        Bz1 = mu0*Hz1
        Bz2 = mu0*Hz2
    return (Hr1+Hr2,Hz1+Hz2,Br1+Br2,Bz1+Bz2)
                    

La fonction suivante permet de dessiner une flèche au milieu d'une ligne.

def fleche(x,y,sens,long,style='k-'):
    n = len(x)//2
    xa = x[n]
    xb = x[n+sens]
    ya = y[n]
    yb = y[n+sens]
    z = (xb-xa)+1j*(yb-ya)
    phi = np.angle(z)
    a = np.pi/2+np.pi/3
    xc = xb+long*np.cos(phi-a)
    yc = yb+long*np.sin(phi-a)
    xd = xb+long*np.cos(phi+a)
    yd = yb+long*np.sin(phi+a)
    plot([xb,xc],[yb,yc],style)
    plot([xb,xd],[yb,yd],style)
                    

Le tracé d'une ligne de champ de l'excitation magnétique se fait en partant d'un point à proximité d'une face, jusqu'à un point à proximité d'une face ou bien jusqu'à sortie d'un rectangle centré sur l'aimant. Il faut remarquer que le choix de la valeur de M et de μ₀ n'a aucune influence sur les lignes de champ.

def ligneH(a,zs,zn,M,ri,zi,sens,rmax,zmax,dmin):
    ds = 0.01*sens
    r = ri
    z = zi
    ligne_z = []
    ligne_r = []
    continuer = True
    while continuer:
        ligne_z.append(z)
        ligne_r.append(r)
        (Hr,Hz,Br,Bz) = champAimant(r,z,a,zs,zn,M)
        H = np.sqrt(Hr**2+Hz**2)
        r += Hr/H*ds                                    
        z += Hz/H*ds       
        if (abs(z-zs) < dmin and abs(r)<a) or (abs(z-zn)< dmin and abs(r)<a) or (abs(r)>rmax) or (abs(z)>zmax) : continuer=False
    return (np.array(ligne_r),np.array(ligne_z))
a=1
M = 1 
zs=-2
zn=2
ds = 0.05
figure(figsize=(8,8))
zmax = 10
rmax = 10
longfleche = 0.2
dmin = 0.1 # distance minimale d'approche des faces
for ri in [0,0.2,0.4,0.6,0.8,1.0]:
    sens = 1
    ligne_r,ligne_z = ligneH(a,zs,zn,M,ri,zn+dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = -1
    ligne_r,ligne_z = ligneH(a,zs,zn,M,ri,zs-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = 1
    ligne_r,ligne_z = ligneH(a,zs,zn,M,ri,zn-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
axis('equal')
xlabel('z')
ylabel('r')
xlim(-zmax,zmax)
ylim(-rmax,rmax)
plot([-zs,zs,zs,-zs,-zs],[a,a,-a,-a,a],'r-')
grid()
title('Lignes de H')
                     

aimant-ligneH.pdf

def ligneB(a,zs,zn,M,ri,zi,sens,rmax,zmax,dmin):
    ds = 0.01*sens
    r = ri
    z = zi
    ligne_z = []
    ligne_r = []
    continuer = True
    while continuer:
        ligne_z.append(z)
        ligne_r.append(r)
        (Hr,Hz,Br,Bz) = champAimant(r,z,a,zs,zn,M)
        B = np.sqrt(Br**2+Bz**2)
        r += Br/B*ds
        z += Bz/B*ds
        if (abs(z-zs) < dmin and abs(r)<a) or (abs(z-zn)< dmin and abs(r)<a) or (abs(r)>rmax) or (abs(z)>zmax) : continuer=False
    return (np.array(ligne_r),np.array(ligne_z))
    

figure(figsize=(8,8))
for ri in [0,0.2,0.4,0.6,0.8,1.0]: 
    sens = 1
    ligne_r,ligne_z = ligneB(a,zs,zn,M,ri,zn+dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = -1 
    ligne_r,ligne_z = ligneB(a,zs,zn,M,ri,zs-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
    sens = -1
    ligne_r,ligne_z = ligneB(a,zs,zn,M,ri,zn-dmin,sens,rmax,zmax,dmin)
    plot(ligne_z,ligne_r,'b-')
    fleche(ligne_z,ligne_r,sens,longfleche,style='b-')
    plot(ligne_z,-ligne_r,'b-')
    fleche(ligne_z,-ligne_r,sens,longfleche,style='b-')
axis('equal')
xlabel('z')
ylabel('r')
xlim(-zmax,zmax)
ylim(-rmax,rmax)
plot([-zs,zs,zs,-zs,-zs],[a,a,-a,-a,a],'r-')
grid()
title('Lignes de B')
                     

aimant-ligneB.pdf

3.c. Force entre deux aimants

On cherche à calculer la force exercée sur un aimant par un autre aimant coaxial. Cette force s'écrit comme une intégrale sur les deux faces de l'aimant :

$$F_z=2\pi M{\int }_0^a{r}_1(B_z(r_1,z_n)-B_z(r_1,z_s))d{r}_1\text{(41)}$$

où $B_z$ est le champ créé par l'autre aimant.

Voyons comment calculer une force sans dimensions. Soit M l'unité dans laquelle les moments magnétiques sont exprimés. Si l'on pose ${\mu }_0=1$, le champ magnétique obtenu est sans dimensions et son unité est :

$$B_0={\mu }_{0}M\text{(42)}$$

Si L est l'unité de longueur, l'unité de la force est :

$$F_0=M{B}_0{L}^2=\frac{{\left({\mu }_{0}M\right)}^2}{{\mu }_0}{L}^2\text{(43)}$$

          
def force2Aimants(a1,zs1,zn1,a2,zs2,zn2,M1,M2,N,mu0=1):
    r1 = np.arange(N)*a1/N
    dr1 = r1[1]-r1[0]
    Fz = 0
    for k in range(N):
        (Hra,Hza,Bra,Bza) = champAimant(r1[k],zn1,a2,zs2,zn2,M2,mu0)
        (Hrb,Hzb,Brb,Bzb) = champAimant(r1[k],zs1,a2,zs2,zn2,M2,mu0)
        Fz += Bza-Bzb
    Fz *= 2*np.pi*M1*dr1*r1[k]
    return Fz
                

a1 = a2 = 5 # en mm
L1 = L2 = 20 # en mm
mu0M = 1 # en teslas
F0 = (mu0M)**2/(4*np.pi*1e-7)*(1e-3)**2
zn1 = 0
zs1 = zn1-L1
zs2 = np.linspace(0.1,20,20) # en mm
zn2 = zs2+L2
Fz = np.zeros(len(zs2))
for k in range(len(zs2)):
    Fz[k] = force2Aimants(a1,zs1,zn1,a2,zs2[k],zn2[k],1,1,100) 
figure()
plot(zs2,Fz)
grid()
xlabel("distance (mm)")
ylabel(r"$F_z/F_0$",fontsize=18) 
                

fig7.pdf

print(F0)
--> 0.7957747154594766